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SUMMARY 

The FORTRAN IV computer program for the design of a two-dimensional supersonic 
nozzle with a sharp-edged throat is presented along with the equations used. The nozzle, 
which has uniform parallel flow at the exit, is designed on the basis of two-dimensional 
isentropic flow of a perfect gas. This program requires as input the exit Mach number 
and the specific-heat ratio. The output yields the nozzle contour for the supersonic por- 
tion. Input and output for a sample case are included. 


INTRODUCTION 

Turbine -driven, hydrogen-fueled, open-cycle auxiliary space power systems have 
recently become of interest. The analysis of such systems in reference 1 indicates the 
possibility of using a supersonic turbine. Proper design methods must be available to 
obtain the highest possible efficiency from a turbine of this type. 

The method of characteristics as applied to the two-dimensional isentropic flow of a 
perfect gas is used for the design of both supersonic nozzles and supersonic rotor blading. 
Such a computerized method for designing supersonic rotor blading for any Mach number 
level and specific-heat ratio is given in reference 2. A method for the design of sharp- 
edged-throat supersonic nozzles that produce uniform parallel flow at the nozzle exit is 
described in reference 3. This type of nozzle, an example of which is shown in refer- 
ence 4, is of minimum length, which is desirable for a compact turbine. References 5 
and 6 present tabulated coordinates for several sharp -edged nozzle designs for a specific - 
heat ratio of 1.4. 

Because of the interest in supersonic turbines for auxiliary space power systems, a 
computer program based on the method of reference 3 was written for the general design 
of sharp-edged supersonic nozzles. The method of characteristics is used in the design, 


and the program is applicable to any Mach number level and specific -heat ratio. This 
report presents the equations used and the FORTRAN IV computer program along with a 
description of the input and output. Sample input and output are included. 


SYMBOLS 


A area 

A* throat area (M = 1) 

k variable index for characteristics of family II 

M Mach number 

M* dimensionless velocity, ratio of velocity to sonic velocity at throat 
mj slope for characteristics of family I 
mjj slope for characteristics of family II 
n variable index for characteristics of family I 

X x-coordinate 

y y -coordinate 

y specific -heat ratio 

ju Mach angle, deg 

/ 

Prandtl -Meyer angle, angle through which flow must turn from Mach 1 to required 
Mach number, deg 

Au increment in Prandtl -Meyer angle, deg 

<p flow angle, angle between velocity direction and x-axis, deg 

I characteristic of first family 

II characteristic of second family 

Subscripts: 
e exit 

k variable index for characteristics of family II 

max maximum 

n variable index for characteristics of family I 
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METHOD OF ANALYSIS 


The design of a two-dimensional supersonic nozzle with a sharp-edged throat describ- 
ed herein was based on establishing parallel uniform flow at the nozzle exit. The method 
of characteristics as applied to the two-dimensional isentropic flow of a perfect gas was 
used for the design in a manner analogous to that described in references 3 and 7. A 
sketch of a typical nozzle designed in this manner is shown in figure 1. 


Subsonic 

portion 



Figure 1. - Sharp- edged- throat supersonic nozzle. 

Since the nozzle is symmetrical about the centerline, only half the nozzle need be de- 
signed, as shown in figure 2. The sharp-edged throat initiates a Prandtl -Meyer expan- 



Figure 2. - Nomenclature and wave diagram for supersonic nozzle with sharp-edged throat. 
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Sion (flow around a corner). The waves or characteristic lines emanatit^ from the sharp 
edge have a negative slope and are termed waves of family II to be consistent with refer- 
ence 7, In the region near the throat, only family II waves exist. These family II waves 
are then reflected at the centerline into waves with positive slopes. These waves are 
termed waves of family I. In the center region of the nozzle, the waves intersect, and 
both families of waves are present. The family I waves extend beyond the region where 
both families exist and intersect the nozzle contour, which is shaped so as to cancel these 
waves. In the exit region of the nozzle, the flow is then parallel and uniform. 

The two physical boundaries for the design procedure are the nozzle contour and the 
centerline. The internal boundaries are the Mach lines or characteristic lines (ref. 7). 
Therefore, each region is bounded by either Mach lines alone, or Mach lines and a 
physical boundary. The equations and calculation procedures needed to design each re- 
gion of the nozzle are given herein. 


Equations 

The nomenclature for the nozzle is given in figure 2. The supersonic portion of the 
nozzle is divided into a finite number of small regions in which the flow properties are 
considered to be constant. Each small region is denoted by two indexing variables k 
and n, where k is a variable index for family II characteristics and n is a variable 
index for family I characteristics. For nozzles with a sharp-edged throat, the contour 
angle at the throat is equal to one-half the exit Prandtl -Meyer angle (ref. 7), where the 
Prandtl -Meyer angle v is defined as the angle through which the flow must turn from 
Mach 1 to the required Mach number M. The relation between v and M is given bv 
(ref. 5) 


M^ - 1^ - arctan - 1 (1) 

The maximum values for k and n are determined by the exit Prandtl -Meyer angle and 
the increment in Prandtl-Meyer angle Au. Therefore, the maximum values of k and 
n are 



k 

2 Au 


( 2 ) 


and 
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’^max 


1 

2 


u 

1 

Av 


( 3 ) 


The flow angle for each small region is 

^ = (k - l)Ai/ (4) 

where k varies from 1 to k^„^. The flow angle is zero near the centerline (k = 1) and 
p^/2 at the wall (k = 

The Prandtl -Meyer angle for each small region ^ is 

^k n ~ 2Ai/(n - 1) + (k - l)At' (5) 


where n varies from 1 to n^„^. The Prandtl -Meyer angle is zero at the throat (k = 1, 
n = 1) and at the nozzle exit (k = 1, n = n_„^). Many of the small regions have the 

0 IlldX 

same Prandtl-Meyer angle, such that ^ n+1' Therefore, the Prandtl-Meyer 

angle used in the program can be designated as Uj, where I = k + 2(n - 1). 

The Mach angle for these increments is 




k, n 


= arcsin 



(6) 


Since v, „ is known, an iterative procedure is required to solve for M, The rela- 
Kj n ^ Kj n 

tion between v and Mf the ratio of velocity to the sonic velocity at the throat, can 
be used. This relation is (ref. 2) 



The bounds on M^ ^ of equation (7) are 1 and V (y + 1)/ (y " !)• Equation (7) can be 

solved for M^ Ijy several numerical methods. The method used in the program is 
n 

Newton's method (ref. 8). The Mach number is then determined from 
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( 8 ) 


M 


k, n 







Thus, the Mach ai^le ^ can be obtained from equation (6). 

The slope of the characteristic lines of family I is tan(jix + <p) and that of family n is 
-tan(ju - <p) (ref. 7). However, since the regions are finite, a better estimate of these 
slopes is the slope based on the average value of the angles of the adjacent regions. 
Therefore, 

m, = tan ^ ^ (i 

V 2 2 ) 

and 




( 10 ) 


The X, y-coordinates of general intersection can now be determined. From analytical 
geometry. 


Vn 


, (yfal.n-l - ■”n^fal,n-l) - (yk-l.n ~ ”l\-l.n) 


m. 


m 


n 


( 11 ) 


and 

yk,n = yk-l,n + >”l(\,n-=^-l,n) (12) 

These equations are general and apply directly to the interior points. The special equa- 
tions for the throat, nozzle contour, and centerline regions are given in the following sec- 
tion along with the procedure. 


Calculation Procedure 

The nozzle design is based on a throat half -height of 1. Therefore, the coordinates 
at the throat are x = 0 and y = 1. The region near the throat where n = 1 requires a 
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special set of coordinate equations. At the point n = 1, k = 1, equations (11) and (12) re- 
duce to 






yi.i 


m 


II 


= 0 


When n = 1 and 2 < k < (k^„^ - 1), equations (11) and (12) reduce to 




_ ^-[yk-i,i-”^iVi, i] 


mj-mjj 


(13) 


(14a) 


y(k, i) = yk-i,i + "^i(^k,i-^k-i, i) 


(14b) 


At the nozzle contour point where n = 1 and k = the contour must be shaped so that 

it cancels the waves. Therefore, 




(15) 


Substituting equation (15) into equation (14a) and then solving equations (14a) and (14b) 
give the nozzle contour point for n = 1. 

For all other values of n, 2 < n < - 1), the following coordinate equations are 

used. At the centerline of the nozzle where k =1, mj = 0. Therefore, equations (11) 
and (12) reduce to 


\n 


yk+l,n-l ~ ”^il^k+l,n-l 


(16) 


m 


n 


yk. n 


= 0 


For the interior points 2 < k < (k^^^ - 1), equations (11) and (12) are used. For the 
nozzle contour points k = k_„„, equations (11) and (12) are used with equation (15) for 

nicix 

mjj. As n increases by 1, k^^^^ is reduced by 1 until k^^,^ = 2. This procedure gives 


max 


all the required nozzle contour points. 
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Increment Size 


The increment in flow turning Av is the choice of the designer. The accuracy re- 
quired determines the increment size used. The exit y-coordinate is known from one- 
dimensional isentropic flow and is equal to the area ratio A^/A*. Reference 5 suggests 
a method of extrapolation to obtain the exit x-coordinate for zero increment size. An 
example of this method is shown in figure 3, where the increment size Av is plotted 
against the exit x, y -coordinates. It is apparent from figure 3, that for an increment 
size of approximately 0. 1°, the errors in the exit x, y-coordinates are 0, 5 percent and 
0. (J2 percent, respectively. For a Ai/ of approximately 0. 5°, the errors in the exit x, 
y-coordinates are 2. 5 percent and 0. 08 percent, respectively. 
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Figure 3. - An example effect 
of increment size on exit 
coordinates. Exit Mach 
number, 4.05; specific- 
heat ratio, 1. 36. 


DESCRIPTION OF INPUT 


A description of the input for the FORTRAN IV computer program is given in this 
section. The input quantities are exit Mach number, increment in flow turniiig, specific- 
heat ratio, and an indicator controlling printout. 
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The input format with sample data is shown in table I. 


TABLE I. - SAMPLE INPUT DATA 
[Numbers In top row are card column numbers.] 


1 6 

7 12 

13 18 

19 20 

ME 

DV 

GAM 

NP 

4.05 

0.10 

1.36 
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The input variables are 

ME nozzle-exit Mach number, M 

’ e 

DV increment in flow turning, Av, deg (this is an initial value; the program will 
adjust this value to make k^„_, an integer) 

m<xX 

GAM specific -heat ratio, y 

NP nozzle coordinate printout indicator (for example, if NP = 1, every coordinate 
is printed out; if NP = 10, every tenth coordinate is printed out; it should also 
be noted that for any value of NP, the last coordinate is always printed out, 
and NP must be an integer, and right adjusted) 


DESCRIPTION OF OUTPUT 

A sample of the output obtained from the program is given in table II. The first line 
is input, except for VE, the Prandtl -Meyer angle based on the input exit Mach number, 
and DELTA V, the calculated increment in flow turning used in the program. The terms 
X and Y are the x, y -coordinates of the nozzle contour. 
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PROGRAM DESCRIPTION 


Main Program 

All the calculations are made in the main program except the calculations for Mach 
number and Mach angle, which are made in Subroutine UA. The program variables for 
the main program are 


I 

k + 2(n - 1) 

KMAX 

integer (eq. (2)), k^^^ 

NC 

index for characteristic of family I, n 

NK 

number of Mach angles 

NMAX 

integer (eq. (3)), n^^^^ 

NMA 

NMAX - 1 

PHI(K) 

flow angle, cp^ 

SLOPE 1 

slope of family II waves 

SLOPE 2 

slope of family I waves 

u(i) 

Mach angle, Uj 

X(K, N) 

x-coordinate 

Y(K, N) 

y-coordinate 

XC(NC) 

x-coordinate of contour point 

YC(NC) 

y-coordinate of contour point 

XT 

x-coordinate at throat 

YT 

y-coordinate at throat 

Subroutine 

UA 


This subroutine calculates the Mach number and Mach angle. Equation (7) is solved 
bv Newton's method for Mf . The Mach number and Mach angle are then determined. 
The calling sequence for UA is as follows: 

CALLUA(I, DV, GAM, U(I)) 
where 

I value of k + 2(n - 1) 

DV increment in flow turning (input), Au 
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GAM specific -heat ratio (input), y 

U(I) Mach angle (output), /Xj 

The program variables for UA are 


DVDMC 

derivative of eq. (7), di^/dM* 

M 

Mach number, M 

MC 

M* 

V 

Prandtl -Meyer angle, v 

VA(MC) 

function statement, where VA is eq, (7) 

VI 

temporary storage location for V 

U 

Mach angle, ju 

X2M 

limit on I., - 0. 01 
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PROGRAM LISTINGS 


SIBFTC NOZZLE DECK 

C TViO DIMENSIONAL SUPERSONIC NOZZLE WITH A SHARP-EDGED THROAT 

OIMENS ION XC( 1500) , YC ( 15 00 ) , PH 1 1 15 00 ) , U ( 3000 )' » 

1X( 1500,2 ) ,Y( 1500,2) 

CCMMQN GAM 
REAL M,MC,ME,MEC 

4 READ(5,3) ME,DV,GAM,NP 

3 FORMAT! 3F6. 3, 12) 

VE=SORT! (GAM+l . 0.) /( GAM-1 . 0 )) *A TAN ( SORT! ( GAM-1 . 0 )/( GAM+ 1.0 » 

1^(ME’«'*2-1,C) n - ATAN(SQRT(ME=«'*2-1.0) ) 

CV= 0V*.C17453 

XT=0.0 

YT=1.0 

KNAX = INT! .5*VE7DV+1. 5) 

CV=VE/(2.0*FLOAT(KMAX-1) ) 

NMAX = KMAX 
WRITE! 6, 251 ) 

251 FORMAT! IHl, 36X, 59HTW0 DIMENSIONAL SUPERSONIC NOZZLE WITH A SHARP-6 
IDGED THROAT) 

VE= VE*57.295a 
CV =DV*57.2958 
WRITE!6,252) VE,ME ,GAM,OV 

252 FORMAT! 1H0,4HVE = , F8 . 3 . lOX , 4HME =, F8 .3 , 10X,7HGAMMA =,F6.3,10X, 
llOHDELTA V = ,F7.5) 

WRITE! 6, 253 ) 

253 FORMAT! 1H0,5X, 1HX,14X,IHY) 

WRITE !6,254) XT,YT 

C\(= DV’f. 017453 
C FLOW ANGLE CALCULATION 

DC 5 K=1,KMAX 

5 PHI!K) = FLOATIK-1 )^0V 

C MACH ANGLE CALCULATION 

NK = 2*KMAX-1 

DO 6 1=1, NK 

6 CALL UA! I , OV ,G AM , U ! I ) ) 

C REGION NEAR THE THROAT 

DO 100 K=1,KMAX 

NC = 1 

N=1 

I = K+2=j‘!NC-l) 

IF !K.NE.l ) GO TO 10 

SLOPEl = -TAN! !U! I ) + U! I+l) )/2.0-!PHI!K) + PHI!K+in/2.0) 

X!K,N)= - l.O/SLOPEl 
Y!K,N) = C.C 
GO TO 100 

10 IF IK.EO.KMAX) GO TO 20 

SLOPEl = -TAN! ! U! 1 )+U! 14-1) ) /2. C-1 PHI IKl + PHI IK+1) )/2.0) 

SL0PE2 = TAN! !U! I )+U! I+l) )/2.0+!PHI !K)+PHI!K-1) )/2.0) 

X!K,N) = ! 1.0-! Y!K-1,N)-SL0PE2*X!K-1,N) ) )/ 

1! SL0PE2-SL0PE1 ) 

Y !K,N)=Y! K-1,N)+SLCPE2*!X! K,N)-X(K-1 ,N) ) 

GO TO 100 
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20 SL0PE1= TAN(PHKK) > 

SLOPE?. = TAN( (U( I ) + Ut I + in/2.0+(PHI (K) + PHI(K- 1 ) )/ 2 . 0 f 
X(K,N) = ( 1.0-(Y(K-l,N)-SL0PE2*XlK-l,N)n/ 

H SLOPE 2- SI. OPEL I 

Y(K,N)=Y(K-l,N)+SLaPE2*(X(K,Ni-X(K-l«N) ) 
X-C(NC)=X(KMAX.N» 

YCINC)= Y(KMAX,N» 

100 CONTINUE 

IF (MOOINC,NP) .EO.O) WRITE(6,2541 XC ( NC ) , YC( NC > 

C REGION DOWNSTREAM OF THE THROAT 
205 N=2 

KMAX=KMAX-1 

NC=NC+1 

DO 200 K=1,KMAX 
I = K+2*(NC-1) 

IF (K.NE.l) GO TO 201 

SLOPEl = -TAN((U(I)+UII+1) )/2.C-<PHI(K) + PHI(K+in/2.0) 

X (K,N)=-( Y (K+1,N-1 l-SLOPEl^XIK+l.N-l) 1 /SLOPE! 

Y(K,NL=0.C 
GO TO 20 C 

201 IFIK.EO.KMAXI GO TO 202 

SLOPEl = -TAN( (U( n + U( 1+ in/2. 0-IPHI (K) + PHI(K+in/2.01 
SL0PE2 = TAN( (U(I )+U( I+U ) /2. 0+ ( PHI ( K) + PH I ( K-1 ) 1/2. 0» 

X (K,N)=( { Y{K+1,N-1)-SL0PEL*X(K+1.N-1) l-I Y(K-1, N1-SL0PE2 
l*X(K-l,Nm/ISL0PE2 - SLOPEl) 

Y(K,N)=Y( K-1,N)+SL0PE2*(X(K,N)-X(K-1,N) ) 

GO TO 200 

202 SL0PE1=TAN{PHI (K) ! 

SL0PE2 = TAN( (IJ{ I ) + UI I+l) ) /2 . 0+ { PH 1 1 K ) + PH I ( K-1 ) 1/2.0) 
XIK.N )=( ( Y(K+1,N-1 )-SLOPEl*X(K+l,N-l) )-| Y( K-1 , N)-SL0PE2 
1*X(K-1,N) ) )/< SLQPE2 - SLOPEl) 
Y(K,N)=Y{K-1,N)+SL0PE2*(X(K,N)-X(K-1,N) ) 
XC(NC1=X(KMA){,N) 

YC(NC)=YIKMAX,N) 

200 CONTINUE 

IF (MOO(NC.NP) .EO.O) WRITEI6. 254) XC ( NC ) , YC ( NC ) 

IF (KMAX.EQ.2) GO TO 250 
DO 204 K=1,KMAX 
X(K. 1)=X(K,2) 

204 Y(K, 1)=Y(K,2) 

GO TO 205 
250 NNA=NMAX-1 

WRITE (6,254) XC ( N MA ) , YC ( NMA ) 

254 FORMATdH , F 1 0 . 5 , F 1 5 . 5 ) 

GO TO 4 
END 
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$IBFTC U 


DECK 


SUBROUTINE UA( I ♦ DV ,GAM ,U » 

COMMON GAM 
REAL MC.M.MCl 

X2M=SORT( (GAM+i.0)/(GAM-1.0l )-.01 

V = FLOAT! I-1I*DV 

IF(V.EO.O.O» GO TO 5 

EXTERNAL VA 

IF (MC.NE.1.0) GO TO 7 

0VDMC=.02 

GO TO A 

7 A= SORT! (GAM+1 .0)/(GAM-1.0n 

B= ( GAM-1 .0)=»MC**2-GAM 
C= ( GAM+1 .0 »/MC=«‘*2-GAM 

DVDMC= A*( GAM-1.0)*MC/ 1 1 . 0-B**21^* I . 5 ) -MC«*(-3 )<■ I GAM + 1 .0 »/ ! 1 .O- 
1C**2')=»=*( .51 

4 MC= (V-VA! MO ) /DVOMC+MC 
IF (MC.LT.X2Ml GO TO 9 
MC=X2M 

WRITE (6f2) MC.NCtK 

2 FORMAT! 1H0,41HLIMIT HAS BEEN REACHED. MC SET = TO X2M,5X«4HMC 

1F8.3,5X»4HNC =, I4t5X,3HK =.141 
9 Vl= VA(MC) 

IF (ABS! V-Vn.LT. .00001) GO TO 8 
GO TO 7 

ft M = SORT!! I2.0/IGAM+1.0) )*MC**2)/(1.C-!!GAM-1.0)/ 

1 ( GAM+1 .0 )*MC**2) ) ) 

GO TO 6 

5 M=1.0 
MC=1.0 

6 U=ARSIN! I .0/M) 

RETURN 

END 


$I8FTC VAN 


DECK 


FUNCTION VA(X) 

COMMON OAM 

A= SORT! (GAM+l.OI/IGAM-l.O) )-1.0 
B= ARSINUGAM-1.0)*X*«2-GAM) 

C=ARSIN( (GAMFl .0 )/X**2-GAM) 

VA=( 3 . 1415926 ) / 4.0 * At. 5 * (( A+ 1 * 0 ) *B+C ) 

RETURN 

END 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, October 18, 1967, 
128-31-02-25-22. 
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*'Tihe aeronautical and space activities of the United States shall he 
conducted so as to contribute ... to the expansion of human knotvl- . 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof/* 

—National Aeronautics and Space Act of 1958 
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